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We study the dynamic process of splitting a con- 
densate by raising a potential barrier in the center of a 
harmonic trap. We use a two- mode model to describe 
the phase coherence between the two halves of the con- 
densate. Furthermore, we explicitly consider the spa- 
tial dependence of the mode funtions, which varies de- 
pending on the potential barrier. This allows to get 
the tunneling coupling between the two wells and the 
on-site energy as a function of the barrier height. More- 
over we can get some insight on the collective modes 
which are excited by raising the barrier. We describe 
the internal and external degrees of freedom by varia- 
tional ansatz. We distinguish the possible regimes as a 
function of the characteristic parameters of the prob- 
lem and identify the adiabaticity conditions. 

PACS numbers: 03.75.-b, 42.50.-p, 32.80.Pj 



I. INTRODUCTION 

Bose-Einstein condensation of an ideal gas 
is typically presented in introductory textbooks 
solely in terms of particle numbers. And quan- 
tum mechanically enhanced number densities were 
the 'smoking gun' observed in the first experiments 
on dilute gas condensates. But many phenomena 
of interacting condensates depend critically on the 
conjugate quantity to particle number, namely the 
quantum mechanical phase One can have 

highly occupied states with or without phase co- 
herence between them, and the presence or absence 
of phase coherence can make a dramatic difference 
in the physical properties of an ultra-cold gas. As 
in the case of a gas held in an optical lattice, which 
can be a superfluid or a Mott insulator depending 
on the strength of the lattice, the onset or loss of 
phase coherence can even be a phase transition || . 

Since the global phase of a condensate is unob- 
servable, the simplest system in which phase co- 
herence can be manifested consists of two states, 
occupied by a large number of bosons. As such 
a system can realistically be approximated by a 
condensate in a double well, it has recently at- 
tracted attention ||7|-p^|. In these works, Joseph- 
son oscillations and the phase coherence between 
two coupled condensate were studied considering 



time independent coupling parameters. In |l6| ] the 
disappearance of the phase coherence between the 
two wells due to a change in time of the tunneling 
coupling has been studied. The coupling param- 
eters and the related phase coherence properties 
depend on the potential barrier between the two 
wells, since in the limit where the barrier is very 
low we have a single condensate and in the limit 
where the barrier is very high we have two com- 
pletely separated condensates. Even if, at least in 
the high barrier limit, it has often been pointed 
out how to relate the coupling parameters (on-site 
energy and tunneling coupling) to the overlap of 
the wavefunctions localised in the two wells, this 
has never really been taken explicitly into account. 
Introducing the spatial degrees of freedom, as we 
did, allows us to relate all that to the potential bar- 
rier in a more than phcnomcnological way and be- 
comes expecially important in the study of the dy- 
namics of the process, because the wavefunctions 
change drastically in time and collective modes are 
excited. 

In this paper we examine this problem and de- 
velop a method which can be extended to the case 
of many wells, in order to encompass the turning 
on of an optical lattice in a condensate, allowing to 
go beyond a Gross-Pitaevkii treatment, as the one 
done for example in JT7| . The physical situation is 
similar to the ones already realized experimentally, 
where a double well potential was created by shin- 
ing a far-off resonant laser beam in the center of 
the magnetic trap (see e.g. Q) or where an array 
of traps was created by on optical standing wave 

Est 




FIG. I. Schematic diagram of the process. 

We are interested in the full dynamics of the 
process: in addition to the phase coherence prop- 
erties, we want to study the excitation of the col- 
lective modes by taking the spatial dependence 
of the condensate wavefunction explicitly into ac- 
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count. The price of including all these effects with- 
out assuming mean field theory is that we must use 
a time-dependent variational approach, choosing 
variational ansatz to describe both the "internal" 
and "external" dynamics (that is, the distribution 
of particles between two motional states treated as 
given, and the evolution of the spatial wave func- 
tions of these states) . This allows us to reduce the 
intractable full problem to a set of coupled differ- 
ential equation for our few variational parameters. 
Although the time-dependent variational approach 
is not guaranteed to be quantitatively accurate, it 
allows qualitatively important processes to be in- 
vestigated, and it has proven surprisingly reliable 
in previous applications to condensate physics . 
Here we use it to derive the coupling between the 
internal and external dynamics, investigate which 
are the conditions under which the two can be de- 
coupled, and identify the typical time scales for 
both. 

In Sec. H we briefly present our two-state model 
and discuss qualitatively the behaviour of t he sys- 
tem that our model can explain. In Sec. III. wc in- 



troduce explicitly the variational ansatz we choose 
to describe the phase dynamics and the collec- 
tive modes. After introducing the time-dependent 
variational principle and deriving the Lagrangian, 
we write the equations of motio n for the varia- 
tional parameters. Then in Sec. Ill E we discuss 



the conditions under which it is an accurate ap- 
proximation to neglect the coupling between inter- 
nal and external degrees of freedom. This decou- 
pling allows us to model both internal and exter- 
nal evolution in a still simpler way. The external 
dynamics is studied in Sec. IV, comparing simple 
analytic models with numerical results. Sec. [v] is 
devoted to the internal dynamics, where the statics 
is studied and analytic estimates obtained. With 
these decoupled studies to guide expectations, the 
full variational equations of motion, with coupled 
internal and external degrees of freedom, will be 
studied numerically in Sec. VI. The results are 



compared with those of the phase model describ- 
ing the relative phase of two weakly coupled super- 
fluids p3,p2). We conclude with a general discus- 
sion of our results and their implications. Two ap- 
pendices define the number difference and relative 
phase operators , derive the phase model Hamil- 
tonian, and show that our variational ansatz ade- 
quately represents the evolution generated by it. 



II. MODEL 

Let us consider the situation in which we have 
N bosons confined in a harmonic trap at zero tem- 
perature. We slowly deform the trap symmetri- 



cally around its center raising a potential barrier, 
until it becomes a double-well potential. We want 
to study the dynamics of the process as well as the 
final state of the bosons. 

We will treat the problem in one dimension. 
This describes the situation in which we have a 
cigar-shaped condensate, and we deform it along 
the most elongated direction. Nevertheless, our 
model could be extendend in a straightforward way 
to 2 and 3 dimensions. This would just lead to 
more complicated equations, but would not affect 
the main results. 



A. Two— mode model 

The second quantization Hamiltonian describing 
the situation we have in mind and in which parti- 
cles interact via a <5-pseudopotential is 

H{t)= j¥{z) (J?- + V(z,t)^(z)dz + 

+ ±g [fr(z)&{z)9{z)#{z)dz, (1) 



where iff is the bosonic field operator and V(z, t) is 
the time-dependent potential which describes the 
deformed trap. Here, g is an effective coupling 
constant which depends both on the s-wave scat- 
tering length and on the atomic distribution in the 
transverse directions. 

We assume that at time t = the atoms are 
in the ground state of this Hamiltonian, and we 
want to determine the state after the potential is 
deformed. This problem cannot be solved even 
numerically (see, for example, p3| and references 
therein). Therefore, we need to consider a simpli- 
fied model which describes the main features of the 
process. If the process is completely adiabatic, the 
final state will be a 'fragmented condensate' with 
half of the particles in each of the potential wells: 
when the two condensates do not interact this is a 
much lower energy state than the one with phase 
coherence, because minimizing fluctuations in the 
relative particle number lowers the energy due to 
intcrparticle repulsion. Such a fragmented conden- 
sate can be regarded as having two entirely inde- 
pendent condensates. If we changed the potential 
very fast, then we would obtain a single conden- 
sate that oscillates in each of the potential wells. 
We would also expect to see collapses and revivals 
in the "condensate phase" pL provided the losses 
are not important |24(] . 

It is clear that the Gross-Pitaevskii Equation 
(GPE) will not give a good description of the split- 
ting process, in principle. This equation describes 
the evolution of a single mode of the condensate 
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(p(z,t), and, therefore, is not valid for fragmented 
condensates. In order to interpolate between the 
Gross-Pitaevskii limit of phase coherence and the 
limit of two independent condensates, one needs to 
consider at least two modes fi^iz, t). Then, one 
can write the state of the atoms as (we assume N 
to be even) 

N/2 

|3(i)>= £ c m (t)\m(t)), (2) 

m=—N/2 

where 

= I I \vac), (3) 

and ai^ are the mode annihilation operators de- 
fined as 

aj(t) = J Vi{z,t)tf{z)dz. (4) 

We have to consider the evolution of the wavefunc- 
tions yi,2 as well as of the coefficients c m , which 
will be coupled and governed by the Hamiltonian 

H{t) = (5) 

ij=l,2 J \ / 

,5 t t 

+ 2 a i«i a / a m 

ijlm— 1,2 

y ^* l {z,t)ip*{z,t)ipi{z,t)^p m {z,t)dz. 

We will refer to the evolution of those wavefunc- 
tions as "external dynamics" and to the one of the 
coefficients as "internal dynamics" . 

In order to find the mode-functions 991,2 we can 
use the variational principle (in the same way as 
one derives the GPE from a Hartree-Fock ansatz). 
However, this also turns out to be very compli- 
cated. A way around that problem is to express 
tpi.2 in terms of some few variational parameters: 
we will use quasi-gaussian functions to describe 
those mode- functions. On the other hand, we 
could also use a variational principle to determine 
the evolution of the coefficients c m . This again 
turns out to be very complicated, so that we will 
also use a Gaussian ansatz for them. Once they are 
known, in order to estimate when the condensate 
is fragmented, we can look at the eigenvalues of 
the single particle density operator corresponding 
to the internal dynamics only, that is, the matrix 

1 ( (a{ai) (a\a 2 ) \ ffi s 
9 N { (atai) (4oa) J ' W 



where {a\a.j) mean the expectation value on the 
state |$). The eigenvalues A± indicate whether we 
have a single condensate (A+ ~ 1, A_ ~ 0), or a 
fragmented one (A+ ~ A_ ~ 1/2). 

As already extensively discussed in the literature 
, [l6| , the two- mode model has a limited validity in 
the case of low barrier, when in principle one is not 
allowed to neglect higher excited modes. However, 
it becomes more and more accurate the higher the 
barrier gets, since in this case the two lower lying 
modes move closer together in energy compared 
to the higher ones. Hence it should allow a good 
description of the splitting process. 

B. Qualitative behavior 

Before presenting the numerical and analytical 
result coming from our analysis, we will briefly dis- 
cuss the qualitative behavior that we expect from 
the model under study. We will show that the 
equations we will derive for the external and inter- 
nal dynamics can be decoupled to a very good ap- 
proximation. That is, we can first solve the equa- 
tions for the external dynamics basically taking 
constant values for the variational parameters de- 
scribing the coefficients c m . Once these equations 
are solved, we can use the corresponding wave func- 
tions yi,2 to calculate the time-dependent coeffi- 
cients for the equations that describe the internal 
variational parameters. In summary, once we have 
solved the equations for the external parameters, 
we are left with a two-mode model with time- 
dependent coefficients which contain all informa- 
tion about the external dynamics: they will define 
the hopping and on-site interaction, whose compe- 
tition determines the phase relation between the 
two modes. 

Regarding the external dynamics, one can see 
two kinds of behaviors depending on the time scale 
r at which the barrier is raised. The important 
time scale with which one has to compare is the 
oscillation period t z in the trapping potentials at 
each time. These periods change by roughly a fac- 
tor of two between the initial harmonic potential 
and the final double well (with our specific choice 
of the trapping potential). Thus, if r 3> t z the pro- 
cess will be adiabatic with respect to the external 
dynamics, which means that ipi^(z,t) will basi- 
cally correspond to the two ground states of the 
right and left wells at the final time. If t <C t z we 
will have collective excitations, in which tpi^(z,t) 
oscillate strongly. In this case, we will have that 
the energy of the condensate E (with respect to the 
ground state energy) has increased, so that it may 
be destroyed. Although we cannot account for the 
disappearance of the condensate within our model, 
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we can estimate when this will happen just by con- 
sidering the fact that under normal circumstances 
thcrmalization will occur, and, therefore, the con- 
densate will be destroyed for E = KbT c > NTilo z , 
where uj z = 2it/t z , Kb is the Boltzman constant, 
T c indicates the critical temperature and E is the 
extra energy in the final state. We find that the 
condensate disappears for r ~ t z . 

Regarding the internal dynamics, there is also 
a time scale in the problem which determines the 
dynamics; this is the revival time r r . Given two 
condensates with an initial well defined relative 
phase, it is well-known that the relative phase first 
disappears (collapse) and then is restored at time 
r r ~ 1/g If t 3> T r the process will be adi- 

abatic with respect to the internal dynamics, the 
phase coherence will be lost during the process and 
therefore we will end up with two independent con- 
densates in each well, with no phase coherence at 
all (note that it makes no sense to talk about col- 
lapses or revivals in this situation). If r <C r r at 
the end of the process we will have two condensates 
with a good phase coherence. In that case, after 
the splitting, collapses and revivals could be ob- 
served provided the particle losses are practically 
absent. 

In summary, we have two important time scales 
in the problem, namely t z and r r . Typically, in 
experiments r r 3> t 2 , so that it will be harder to 
be adiabatic with respect to the internal dynam- 
ics than to the external one. On the other hand, 
since t t is very long in practise, it will be hard to 
achieve r 3> r r within the validity of our model, in 
which particle losses and other imperfections are 
not included. 

Finally we notice that the external and internal 
dynamics depend in a very different way on the pa- 
rameters g and TV. Very similar to what happens 
in the GPE, the equations of motion for the param- 
eters describing the external dynamics depend al- 
most only on the product gN. On the contrary, the 
on-site energy splitting scales like g and the tun- 
neling coupling like N giving rise to very different 
internal dynamics. For increasing N and decreas- 
ing g the relative phase becomes better defined and 
the time r r required to destroy the phase coherence 
is longer. In particular in the limit N — > oo, unless 
the tunneling coupling exactly vanishes, the GPE 
case of phase coherent condensate is recovered. 



III. VARIATIONAL ANSATZ 

In this Section we introduce a variational ansatz 
to describe the internal and external dynamics. 
To describe the ground state of the system, which 
we will call equivalently static or equilibrium solu- 



tion, two parameters are sufficient: zq which corre- 
sponds roughly to the center of the mode functions, 
and p which is related to the width of the number 
distribution. To allow dynamic evolution we have 
to introduce also the variables 07 and xi, which 
will vanish in the static case. For simplicity in 
the following we present the variational ansatz for 
the symmetric case, describing a spatially symmet- 
ric double-well potential and a symmetric atomic 
distribution among the two modes. Starting from 
symmetric initial conditions (no unbalance in the 
population of the two modes and symmetric mode 
functions (fi(z) = (p2(—z)), the symmetry will be 
preserved under the time evolution. Nevertheless, 
it is possible to generalize our ansatz to describe 
asymmetric situations. We will briefly discuss the 
corresponding results in the conclusions. 



A. Variational ansatz for the mode functions 

For a single condensate in a harmonic trap, a 
Gaussian ansatz for the wavefunction has proved 
to be very useful and able to predict the excita- 
tion frequencies within a very high precision po| . 
In our case we make a similar choice. For a very 
high barrier we expect to find two separate conden- 
sates for each of which a Gaussian ansatz should 
be good. We then define the two functions for the 
left and right side: 



{-) 6XP 



(7) 



exp 



/ pi 2 
—1 — z 
V 2 



In the case of low barrier (small zq), these two 
functions have to be orthogonalized to satisfy the 
orthonormality property required for the two mode 
functions. Hence we define 



L P±{ Z ) = — - 



(tp R ±ip L ), (8a) 



(8b) 



which are different from tpRx because of the two 
different normalization constants for ip + and ip_ . 

The variational parameters describing the mode 
functions 1^1,2 are zq and 07. Physically, the ex- 
citation modes that these parameters can describe 
depend on whether the two Gaussians significantly 
overlap in space or not. If they do, as in the case 
of low barrier, then the excitation mode (changes 
in zq) corresponds to a breathing mode. If they 
do not overlap, as for high barrier, then it corre- 
sponds to an oscillation mode in each of the poten- 
tial wells. Allowing also or to vary (and adding 
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an extra variable for the dynamics) one can de- 
scribe more excitation modes. Instead we fix it to 
a constant value, since the curvature of the poten- 
tial wells will be chosen to be always of the same 
order of magnitude. The value of <jr is not obvi- 
ous since the overlap integrals depend strongly on 
it. For instance the hopping terms can be over- 
estimated due to the long Gaussian tales. To fix 
<tr, we identified a range of values for which the 
dynamic behaviour of the system was qualitatively 
the same and choose one of the values within this 
interval, which turned out to be lower than the one 
corresponding to the static solution. 

Since the mode wavefunctions are linear combi- 
nations of Gaussians, if we choose a trapping po- 
tential of the form 

V(z) = X -Mu? T z 2 + V (t) exp[-^/a(t)], (9) 

the integrals in Eq.(|27j) can be performed analyti- 
cally. 

In the following, we will use dimensionless units: 
dfto = y/h/MwT = 1, hu>T — 1 and h = 1. So, 
all lengths will be measured in units of harmonic 
oscillator length, all energy in units of the trap 
frequency and all times in units of uj^, 1 . 



B. Variational ansatz for the coefficients 

We also take for the coefficients c m (t) a Gaussian 
distribution centered at m = 0, 



= Af(p) exp 



1 

- — h ixi I mr 
4p 



(10) 



where Af(p) is a normalization constant that de- 
pends on p only. The variational parameters are p 
and its conjugate one x/. The parameter p is di- 
rectly related to the width of the distribution |c m | 2 , 
whereas Xi contributes to the width of the Fourier 
transform of such a distribution (i.e., to the width 
of the phase distribution, see App. |Aj). 



C. Time-dependent variational principle 

We study the dynamics using the time depen- 
dent variational principle. To derive the equations 
of motion one starts by writing the action S 



S 



dt- 



($1$) 



2i 



(ii) 



where H and |$) have been defined in (||) and (|2|), 
respectively. In evaluating the term ( < I , | < I > ), one 
should remember that the state |<fr) depends on 



time both through the coefficients and the mode 
functions contained in \m): 



($1$) = 



C rn' 



c m (m m 



(12) 



The Lagrangian which follows takes the form 
1 



L 



2 1 



m 

+ £< a i a j) J <P*<Pjdz-h.c 



(13) 
U. 

(14) 



where 7i= |$) is given by 

n = 

= ( a « a j) / Pi(*>*) (tTm +V(z,t)\ ipj(z,t)dz 
y=i,2 J v ' 

+ \ g Y {a\ a ]aia m ) 



ijlrn — 1,2 



V* («i ( z > *)<M^ t)ip m (z, t)dz. 



From the Lagrangian (|lg|), in general one derives 
the equations of motion, carrying out the variation 
with respect to the discrete variables c m and the 
fields <pi. In our case, all the integrals and expec- 
tation values are functions of the variational pa- 
rameters which can be calculate analytically. The 
only important overlap integrals containing time 
derivatives are 



i ( 1 



2 \<TR 1 

ip 12 ip2,\dz = - ttzzzz^i- 



-2<trzJ 



21 



-2ffR2 



If one defines 

N lt> = (a\a 2 ) + {a\ai) 

and 



&i, (15a) 
(15b) 

(16) 



l(p,Xi, Zq) = 

N f 1 

: "2 w + r 



(17) 



Ntf, z%e- 2aRZ o 
2 1 - e -2°a4 



e -2a R z£ 

the Lagrangian becomes 

L = pxi + l&i — H, (18) 
and the corresponding equations of motion are 
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/ an \ 

/ dp \ 

dH 




These equations describe the internal and external 
coupled dynamics of the splitting of the conden- 
sate. In what follows, we have solved them numer- 
ically in different regimes. Before presenting the 
results we will show that one can decouple the evo- 
lution of the external and internal variables, which 
helps to understand the dynamics. For that, we 
will introduce in the next subsection some analyt- 
ical approximations to derive explicit formulae for 
the quantities TL and X by replacing the discrete 
distribution c m by a continuous one, and treating 
the index m as a continuous variable running from 
— oo to co. 



D. Continuous limit 

In order to calculate Ti. and X we have to eval- 
uate expectation values of the form and 

(alajaia m ) . We can do that if we replace the sums 
in m by integrals extended from — oo to oo. When 
this replacement is valid, we can even calculate the 



width of the number distribution |c„ 



a m , as well 



as the one corresponding to the phase distribution, 
(see App.[A|). We obtain 



VP; 



ff0 = V^ +4pa:? ' 

On the other hand, we have 



( a i 2 2 a i,2) 



N 
iV 2 



P, 



(o{a 2 ) = (ajai)* = 
N ( *f 

= Y exp -y 



2p_ 

N 2 



N 2 



{a\a\aia 2 ) 

(afa 2 ) = 
N 2 



iV 2 
■ exp (-2a? 



P- 



(20a) 
(20b) 

(21a) 

(21b) 
(21c) 



(21d) 
(21e) 



4p 
TV 2 



G^xj 1 



N 2 



Furthermore, in this limit we can determine the 
eigenvalues of the single particle density operator 
corresponding to the internal degrees of freedom, 
obtaining 



A-, 



1 ± exp (~a 2 /2) 



(22) 



Notice that A+ — > 1 for p ~ N/4, xj ~ 0, i.e. ~ 
0, which corresponds to the Gross-Pitaevskii limit ; 
instead A± — > 1/2 for large a^, giving a signature 
of the fragmentation of the condensate. 

One can easily determine the limits of validity of 
this continuous approximation. On the one hand, 
the distribution in c m has to be sufficiently broad, 
which implies p 1. On the other hand, xi has 
to be such that ^ n. For our numerical simula- 
tions, we corrected the expressions in Eqs.(|T]) to 
make them valid Vp and Vx/ : for p > 1 we included 
the periodicity in xi and for p < 1 we performed 
the exact sum over m, considering only the few 
number states which are populated. 



E. Decoupling between external and internal 
dynamics 

The coupling among the Zq dynamics and the 
p dynamics appears in the off-diagonal blocks in 
Eq.(|l9|) and in the dependence of H on all varia- 
tional parameters. In the following we will analyze 
under which condition it is possible to decouple the 
internal (p,xi) and external (zo,07) dynamics, so 
that one can study them independently one from 
the other. 



1. External dynamics 

One can rewrite the equation of motion for zq is 
a more explicit form as 



dX 

dz 



zo 



dH 1 zgexp(-cr fl zg) 
dui 2 1 — exp(— 2ctrz 2 ) 



N*. (23) 



Let us first see under which circumstances it is 
possible to neglect the off-diagonal blocks: 

(i) in the low barrier limit, i.e. exp(— <jrZq) ~ 1, 
the off-diagonal blocks can be neglected if N^, plays 
no role: N^ evolves at the frequency which gov- 
erns the internal dynamics lo v (remember defini- 
tion ([l6])), while the external dynamics is gov- 
erned by the frequency lo z , usually of the order 
of the trapping frequency lot- If u> p ^> u) z , then 
the time average of N^ vanishes and has no ef- 
fect on the evolution of zq. On the other hand 
Lo p ~ gNUi exp(— ctr2: 2 /2) < lo z only for a chemi- 
cal potential \i ~ gNU\ < lot- This happens either 
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for very few particles, a case in which our model 
does not hold, since N ^S> 1 is a fundamental as- 
sumption in our model, or for very weakly interact- 
ing particles, in which N^ ~ 0, because the phase 
coherence is very difficult to be destroyed. 

(ii) in the high barrier limit, normally ia v -C lo z 
(see Fig. |^, small p limit) and N^ might vary 
abruptly, since this is when the phase coher- 
ence is supposed to desappear. Anyway, since 
exp(— (JrZq) ~ 0, the off-diagonal block can be ne- 
glected. 

Now let us analyze the p-xi dependence in H 
(see Eqs. (|l4|,pl]) ) . In the on-site terms this de- 



Finally in Sec. VI, we show the results for the in- 



pendence is of order p/N 2 ,xi/N 2 -C 1 an d can 



be therefore safely neglected (Eqs. ( |21a| , 21q )) . In 
the hopping terms, it scales like exp(— cr'^/a) (a = 
2,1/2, see Eqs.(21c,|21c)): it is strong only when 
the hopping terms are already small and negligible 
in comparison with the on-site terms, so it can also 
be neglected. 

After these considerations, we conclude that the 
zo-dynamics is, in a good approximation and 
in reasonable regimes, independent of the p- 
dynamics. This is confirmed by the numerical sim- 
ulations. 



2. Internal dynamics 

The off-diagonal blocks can be neglected in the 
p-dynamics if dTi/dzo ~ or exp(— ctrz 2 ) ~ 0. If 
the barrier is raised starting from a condensate in 
the ground state, zq evolves almost adiabatically 
in the low barrier limit (dTi/dzo ~ 0); when it 
reaches the high barrier limit, then exp(— orz 2 ) ~ 
0. Therefore, the off-diagonal terms can be ne- 
glected during the whole process. 

Instead, the zq — 07 dependence in Ti is strong. 
One can solve the complete coupled dynamics or 
substitute in "H the adiabatic solution for zq and 
compare the results. We will show some examples 
in the following and see that the difference is small. 



ternal dynamics, which are also useful as a check 
of the decoupling assumption. 

The fact that the internal and external dynamics 
decouple under the condition discussed above, al- 
lows us to estimate the excitation of the collective 
modes using a very simple model. We point out 
that in this case, similarly to the Gross-Pitaesvkii 
equation, the external dynamics is governed by the 
product gN and not by the two quantities sepa- 
rately. 

We choose to raise the potential barrier with the 
following time dependence 



Vo(t) 



V 0fi 



Uanh - + 1 



(24) 



From the static solution, obtained by minimizing 
Hip, Xi, Zq, <jt, Vq) with respect to all variational 
parameters simultaneously at fixed Vq, we know 
the equilibrium position of z at any value of the 
barrier. The zq dynamics can be approximated 
very well by the dynamics in a harmonic potential 
whose center moves from zoi n to zq fi n following the 
adiabatic solution and whose frequency changes 
corresponding to the frequency of the small oscil- 
lations. 

To get analytic estimations, we model this dy- 
namics by fixing the frequency and shifting the 
center of the potential along a trajectory z c (t) for 
which an analytic solution exists. When the center 
follows an hyperbolic tangent trajectory with time 
constant r, the semi-amplitude of the oscillation is 
given by 



Sz = 
(zofi 



1" /TTLUT\ /TTLOT\ 

zoin) — — csch ( — — J scch ) 



(25) 



4 4 

Otherwise, if the center moves according to a lin- 
ear ramp with time constant r, the average semi- 
amplitude of the oscillation is 



Sz = V2 



ZQfin — ZQi. 



(26) 



IV. EXTERNAL DYNAMICS: 
EXCITATIONS 

In the previous section, we have written the 
equations of motion describing the full coupled dy- 
namics and demonstrated that if one splits a con- 
densate starting from a ground state configuration, 
it is possible to decouple the internal and the ex- 
ternal dynamics. In this section, we will use this 
result and study the external dynamics decoupling 
it from the internal one. Making use of the exter- 
nal static solution, in Sec. [v| we will discuss the 
static solution for the internal degrees of freedom. 
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FIG. 2. (a) Static solution for 20 (full line) and po- 
sition for the center of the harmonic potential for a 
raising time scale t: hyperbolic tangent shift (dashed 
line) and linear shift (dotted line); time derivative of 
zq (t) (inset), (b) Semi-amplitude of the final oscilla- 
tions (left axis) and final excitation energy per particle 
(right axis): numerical results (•), results for the tanh 
(□) and linear shift (o) for lo = 2.82 and 1.54. 
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We compare in Fig.g these expressions with the 
numerical results. For fast raising of the barrier the 
hyperbolic tangent shift gives a good estimate. For 
slower raising, instead, the amplitude of the oscil- 
lations is largely underestimated. In this case, the 
linear shift of the center is useful to give an upper 
bound. We have checked numerically (by chang- 
ing the frequency of the harmonic potential in time 
according to the adiabatic solution) that the small 
discrepancy between the actual shift of the center 
and the hyperbolic tangent dependence is enough 
to produce such a big change in the amplitude of 
the final oscillations. The reason for this might be 
better understood comparing the time derivative 
of z c (t) in the two cases (see inset in Fig.||(a)). 
Our estimations are qualitative, but allow us to 
set some lower and upper bounds and deduce an 
adiabaticity condition r <C 1/lo z for the external 
degrees of freedom. In the same way we can es- 
timate the extra energy per particle due to the 
excitation of collective modes. Using the relation 
KbTc ~ NTwj z (the trapping frequency at the end 
of the process is oj z and not wt), we can estimate 
which is the fastest time scale which does not de- 
stroy the condensate. As expected, it is r ~ 1/uj z . 



V. INTERNAL DYNAMICS: STATIC 
SOLUTION 

In this section, we will study the complete split- 
ting process, starting from a condensate trapped in 
a harmonic potential and ending with a potential 
barrier with height Vq- Apart from presenting the 
numerical results obtained by solving Eqs . jl9|) , we 
will also introduce a simple two-mode model that 
allows a deeper understanding of the results ob- 
tained by the variational ansatz. 



A. Two— mode model 

Given the fact that the external dynamics is ba- 
sically independent of the internal one, we can 
derive a simple model that accounts for most of 
the effects related to the internal dynamics. The 
Hamiltonian in Eq.(|5|) depends on the following 
overlap integrals 

f P 2 

K H = J Vi{ z )2M<Pj(z)dz, (27a) 

Vy = j ^(z)V(z)<pj(z)dz, (27b) 

Ux = J \<p t (z)\ 4 dz, (27c) 

U 2 = f \<p i (z)\ 2 \<p J {z)\ 2 dz = (27d) 



r,= I \cp. l (z)\ 2 i P *(z)cp J (z)dz, (27e) 



where i, j =1,2 p.%. We use 
a\a\aia,2 + a\a\a 2 a 2 = 



(28) 



a\ai - 1 + a\a 2 ) a\a 2 « Na\a 2 , 



to define two effective single-particle hopping 
terms J± 2 = —K X2 - V\ 2 - gNU 3 = —K21 - 
V 2 i — gNU 3 . The static solution for the exter- 
nal dynamics is known from the minimization of 
li.{p, xj, Zo, 07, Vo)- P lug ging the corresponding so- 
lution zo(Vo) in Eqs. (071), we get Hamiltonian pa- 
rameters depending only on the barrier height Vo 
(in particular J\ 2 = J 2 \ = J). Neglecting constant 
terms, we write the simplified Hamiltonian 



|2 2 1 t2 2 

a[ a 1 + a 2 a 2 



- J 



a\a 2 + a 2 a\ 



2gU 2 a\a\a 1 a 2 + ^gU 2 



|2 2 1 t2 2 

a[ a 2 + a 2 a l 



(29) 



As it is well-known (see App. |a|), under certain 
conditions we can replace this model Hamiltonian 
by a phase model of the form 



H = -g (Ux - 2U 2 )^-JN cos < 
TV 2 

+9— U 2 cos2cj), 



(30) 



where (f> represents the relative phase between the 
two modes. The overlap integral U 2 may be non- 
negligible at the beginning of the process when the 
two mode functions overlap in a sensible way. In- 
stead at the end, when the two condensates are 
almost spatially separated it is very small. In this 
case U 2 — > and one recoveres the Josephson's 
Hamiltonian ]25j . The ground state of such Hamil- 
tonian is a localised wavefunction for JN ^S> gU\ 
(corresponding to a broad number distribution) 
and a delocalised one for JN <C gU\ (correspond- 
ing to a narrow number distribution). 



B. Static solution and check of the Gaussian 
ansatz 

For this simplified two-mode model, we write 
analytic approximated espressions for the static 
solution and check the validity of the Gaussian 
ansatz. As explained above, we let the parame- 
ters J, Ui, U 2 depend on the barrier height Vo, 
according to the static solution. The expecta- 
tion value of the Hamiltonian can be now written 
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as a function of the internal degrees of freedom 
only, H(p,Xi,Vo), and allows to study the coher- 
ence properties of the ground state at the differ- 
ent stages of the splitting process: for increasing 
barrier height, we expect the fluctuations in the 
number distribution to become smaller. 

We solved the static problem numerically, find- 
ing the minimum of H. with respect to p and Xi 
for fixed Vq. Moreover, in the limits of large p 
(exp(— er^/2) ~ 1) and small p (p — > 0) it is pos- 
sible to get analytic estimation for the value of p 
at equilibrium and the frequency uj p of the small 
oscillations as a function of the overlap integrals 
pijjfl. In the large p limit, one finds 



2( J - gNU 2 ) 



4 V 2 J - gNU 2 + gN{Ux - 2U 2 ) ' 



(31a) 

Lu p = (31b) 
2^2(J - gNU 2 )[2J - gNU 2 + gN^ - 2U 2 )]; 

in the small-p limit (JN <C gU), where the contin- 
uum approximation for m is not valid, p s and io v 
can be calculated with perturbation theory, consid- 
ering the Hamiltonian in Eq.(29) with J = U 2 = 
as unperturbed Hamiltonian and the number state 
| m = 0) as unperturbed ground state. Then we 
get 



4 log 



2gUi 
JN 



gUi. 



(32a) 
(32b) 



In Fig. ^ we plot the two analytic solutions for 
p s , comparing them with the numeric solution for 
the minimum of TC using our ansatz and the cor- 
responding value of p obtained by minimizing the 
exact Hamiltonian for N = 200 and g — 5. The 
agreement between the variational solution and ex- 
act one is excellent and the analytic expressions 
interpolate correctly in the limits where they are 
expected to work. 



Given the analytic espressions in Eqs.([U 32), we 
observe that the oscillation frequency for large p 
coincides with the breathing frequency in the har- 
monic potential approximating the cosine poten- 
tial, and the oscillation frequency for small p cor- 
responds to the revival ti me wh en the cosine poten- 
tial is negligible (see Sec. VI A ) . Aware of the fact 
that it is not possible to define a transition point 
between the two regimes, being it a smooth tran- 
sition, we still calculate the value of p at which the 
two frequencies coincide, get p = 0.125 and claim 
that the phase relation between the two conden- 
sates is smeared out for p < 0.125. 



VI. RESULTS: DIFFERENT REGIMES 

In this section, we show some numerical results 
obtained by integrating the equations of motion for 
the variational parameters. In those results the full 
coupled dynamics of the process were considered. 
It is possible to compare with the evolution of the 
phase distribution governed by the Hamiltonian 
( |30| ) with time-dependent coefficients. Moreover, 
it is possible to get analytic estimates concerning 
the typical time scales of the process. 

In the three cases presented below we will fix 
the product gN in order to have similar external 
dynamics and better isolate the effect of differ- 
ent g and N on the phase properties of the sys- 
tem. At the end of the process, when the bar- 
rier has reached its final value, depending on N 
and g, one can have gUi 3> JN or JN ^> gU\ 
(we assume that U 2 is then negligible). The case 
gU\ ^> JN corresponds to the situation where the 
splitting process is completed and one expects a 
ground state with no well-defined relative phase. 
Instead in the case JN ^> gU\ the ground state 
is still characterized by a localized phase distribu- 
tion, and even if the two condensates are almost 
spatially separated they cannot be considered as 
independent. In this sense, the splitting is not 
complete. 
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FIG. 3. Static variational solution for p as a func- 
tion of a (full line), analytic solutions in the large and 
small p limit (dashed-dotted line), solution of the mini- 
mization of the exact Hamiltonian (circles). All results 
are for N = 200 and g = 5:a parametrize the barrier 
(see Eq.(|) with V = ab §§). 



A. Complete splitting 

We first analyse the case where in the final stage 
of the process, one has gU\ 3> JN. Since in this 
subsection and in the following we fix the product 
gN, this case is obtained for relatively small N and 
large g. Depending on the time scale of the process, 
it is possible to observe collapses and revivals or to 
reach a final fragmented condensate characterized 
by very small number fluctuations. 

We assume that in the final stage of the process 
it is possible to neglect JN and N 2 U 2 , since they 
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depend on the overlap of the two mode functions. 
Then the eingenstates are the number states 



H\m) = gUim 1 



(33) 



and the time evolution of the final state corre- 
sponds to 



!*(*)) « 

^exp 



77? \ 

.—jeKpH^mt) 



(34) 



In the variational ansatz formulation, one can plot 
the constant energy trajectories for the final bar- 
rier height (see Fig. |](a)). Then, the time evo- 
lution corresponds to p and xj following one of 
these trajectories: p keeps almost constant and xj 
evolves unbounded increasing linearly with time 
with "velocity" g\J\ (see Eq. (|32b[) ), which exactly 
reproduces the phase in Eq.(t 
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FIG. 4. Contour plot of H for a = 120: 



TV = 2 x 1(F, g = 0.5, (b) TV 
and (c) N = 2 x 10 5 , g = 0.005 . 
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The main features of this time evolution are the 
following: the width of the number distribution is 
constant (determined by Pfi n ) and do not evolve in 
time; instead the phase distribution collapses and 
revives: an initial distribution peaked around zero 
smears out, revives around (j> = it and so on. Defin- 
ing the collapse time as the time r c when cr^ ~ n 
and the revival time as the time r r when the orig- 
inal phase distribution is recovered shifted by 7r, 
one gets EM. 



(35a) 
(35b) 



Apf in gUi 
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The collapse time is governed by Pfi n , which de- 
pends in general on the barrier raising process and 
on the parameters g and TV. 

The final width of the number distribution ob- 
tained after the raising of the barrier is completed 



depends on the time scale of the process. This is 
just given by the value Pfi n at which the number 
fluctuations are frozen. To evaluated it, we claim 
that as long as uj p > 2t:/t the number fluctuations 



follow the static solution, and when 



2w/r, 



they are frozen out to a final value Pfi n - Setting 



2tt/t in Eq.(??) and substituting in Eq.(31a 



one gets [E9| 



TV 



Pfi 



1 



2tt 



4 2 [2 J - gNU 2 + gN{Ux - 2U 2 ) 
1 2tt 



8gU! t 



(36) 



Of course Eq. (pq) is not valid for r — > oo. For 
t > 2tt/ gU\ = 2r r , we are in the adiabatic regime: 
uj p > 2tt/t during all the process, and we expect 
to reach a completely delocalised relative phase. 

We check this numerically and plot the results 
in Fig. H for several different splitting processes 
(TV = 2 x 10 2 and N = 2 x 10 3 ). We actually find 
that for time scales r > 2r r the final p values lie 
in the region p < 0.125. For faster time scales, 
the estimation in Eq.(j3^) was shown to be very 
good when the external degrees of freedom evolve 
adiabatically. Otherwise discrepancies can be ob- 
served (Fig. H). It is not straighforward to estimate 
such discrepancies, since they depend on the exact 
dynamics and can be either positive or negative. 
Anyway, they are not striking and do not change 
from a situation in which the final relative phase 
is well defined to the opposite one. 




r, 2tt/uj. 



p 



FIG. 5. Case (i): TV = 2 x 10 2 , g = 5 and 
TV = 2xl0 3 ,<? = 0.5. Final p value against r: com- 
plete dynamics (•) and adiabatic external dynamics 
(>); analytic solution for pfm against 2h/uj p in the 
large and small p regimes (full lines) ; the shaded area is 
for p < 0.125 and represent the situation of completely 
delocalised relative phase; Case (ii): TV = 2 x 10 s , 
g — 0.005; analytic solution for p/i„ (full lines) and nu- 
merical value for pu (o); Static solutions for a = 120 
(o) H). 



VI RESULTS: DIFFERENT REGIMES 
1. Collapses and revivals of the phase 

Now we consider two of these splitting processes 
more in detail and compare quantitatively with the 
evolution of the phase distribution following the 
Hamiltonian in Eq. (|30|) , where J, U% and U2 vary 
in time with the barrier height according to the 
instantaneous static solution. 

We have already mentioned that collapse and 
revival of the phase are predicted for two conden- 
sates with an initially good phase relation when 
the final tunneling coupling is negligible. Let us 
consider the case N = 2 x 10 3 , g = 0.5 and r = 4. 
In Fig. |^ we show the one-atom density matrix 
eigenvalues and the indetermination in the phase 
distribution o^. The agreement between the re- 
sults of the two simulations is perfect [|o| and our 
analytic estimations are also confirmed: we expect 
t c ~ 8 and t t ~ 50, which agree very well with the 
numerical results shown in Fig. ||. 
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FIG. 6. Eigenvalues X± (a) and o^/n (b). Numer- 
ical solutions for the variational ansatz (full line) and 
for the phase model (dotted line); N = 2 x 10 3 , g = 0.5 
and t = 4. 

The actual possibility of observing the revivals 
of the phase in an experiment is something outside 
our model. If they are actually destroyed by par- 
ticles losses p4|] , one is left with two condensate 
with no phase relation, but higher number flucta- 
tions than they would have in the ground state. 

2. Final fragmented condensate 

Another way to cut the initial condensate into 
two independent ones, is to raise the barrier much 
slower, so that the phase coherence is lost adia- 
batically all along the process. So, we now choose 
again iV = 2 x 10 3 , <? = 0.5 but a longer time scale 
r = 200. 

The agreement between the results of the vari- 
ational ansatz and the phase model is very good 
also in this case. The final state is characterised 
by much smaller number fluctations with respect to 
the previous case (see Fig. 0(a)) and by a complete 
delocalised relative phase as shown in Fig. R(b). 
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FIG. 7. Eigenvalues A± (a) and a^/n (b). Numer- 
ical solutions for the variational ansatz (full line) and 
for the phase model (dotted line); N = 2 x 10 3 , g = 0.5 
and r = 200. 

From our analytic estimations, we actually ex- 
pect to reach the static solution for r > 2r r ~ 100. 
Instead as it can be seen in Fig. | for N = 
2 x 10 3 ,r = 200, this does not happen. Dealing 
with such small values of Pfi n it is very easy to 
get at the end very small excitations, which can 
be due both to the external degrees of freedom or 
to some excitations already present in the initial 
conditions. The important feature is that the final 
relative phase is anyway completely delocalised. 

B. Incomplete splitting 

We now analyse the case where in the final stage 
of the process, one has JN ^> gU\. This case is 
obtained for large N and small g: in the limit of 
large N, even if J might be very small, it can hap- 
pen that JN > gU\ and the cosine potential in 
the phase representation (Eq.(pO|)) is not negligi- 
ble. This can be seen as a process in which the 
barrier is raised up to a level at which the splitting 
is not really completed. 

In the case in which the cosine potential at the 
end of the raising process is still deep, so that the 
lowest levels can be approximated with harmonic 
oscillator levels spaced by y / 2JNgUi, the time evo- 
lution follows 

oc ^]c„exp (-iny/2JNgUxt \ \n), (37) 

n 

where \n) are the harmonic oscillator eigenstates 
and where the coefficients c n depend on the exact 
dynamics of the raising process. In particular, for 
symmetric initial conditions, the phase distribu- 
tion is always symmetric and only the even eigen- 
states are populated. Then, the phase distribution 
breathes with a frequency 2\/2 J N gU\, remaining 
always centered at <f> = 0. Moreover we notice that 
in such a case, the width in the number distribu- 
tion is not expected to be constant. In the varia- 
tional ansatz treatment, we have to look again at 
the orbits in the phase space. The contour plot 



VII DISCUSSION AND CONCLUSIONS 



12 



in Fig. ^(b) is just to show how the orbit modify 
from the limit of small N to the limit of large N. 
So let us consider Fig. |I](c). The orbits around 
the minimum of Ti. represent a time evolution in 
which both the width of the number and phase 
distribution change in time. The frequency of the 
small oscillations around the equilibrium position 
can be calculted analytically for exp(— er^/2) ~ 1 
(see Eq.(??) for U2 = 0) and coincides with the 
breathing frequency of the phase distribution in 
the case of the superposition of harmonic oscillator 
eigenstates if gNU\ 3> J. If this condition is not 
satisifed, one is not in the weak coupling regime 
and the phase model is not valid. 

The orbits in the p-xi space are characterised 
by very large oscillations in p. Hence we cannot 
talk of frozen number fluctations. Nevertheless, 
with arguments similar to the ones used before, 
we can try to identify the orbit which describes 
the dynamics at the end of the process. We es- 
timate the maximum p value in such an orbit to 
be pm ~ Pfin- The agreement with the numerical 
solution can be checked in Fig. ^ and it is within 
a factor of 2. The adiabaticity condition in this 
case consists in requiring that the final state is 
superfluid, as the static solution would be. This 
means that the minimum value of p correspond- 
ing to the same orbit as pm must be such that 
the phase coherence is still good. We found that 
a final phase coherence corresponding to a mini- 
mum of A + = 1 — (3 (with (3 <C 1/2) is reached in 
process with typical time scales Tp — 2ir/8JN(3. 
Note that this condition is weaker than requiring 
Pfin — Ps- In fact rp < t s = 2ir/ \/8N JgUi, since 
(3 < \J 'gUi/8JN , as it would correspond to the 
static solution. This means that we still allow even 
big oscillation of p around the equilibrium value, 
as long as they do not destroy the phase coherence. 
This corresponds to a breathing of the phase dis- 
tribution which never becomes completely smeared 
out. Moreover, while r s depends only on the prod- 
uct gN, the adiabatic time scale rp scales like 1/N, 
getting easier and easier to be met for large N. 



1. Final superfluid phase 

As done before, we take now a particular case 
and check directly the results with the one ob- 
tained in the phase model. We choose N = 2 x 10 5 , 
g = 0.005 and t = 50 and find that the eigenval- 
ues A± oscillate with frequency ui p (see Fig. ||(a)): 
A + is always close to 1 and A_ close to 0. This 
corresponds to the breathing of the phase distri- 
bution in the non negligible cosine potential or in 
the variational ansatz treatment to one of the or- 
bits in Fig. |J(c). No complete smearing out of the 



phase is observed (see Fig. ||(b)). Those results are 
confirmed up to a good level by the phase model 
(see Fig. ||(a,b)), even if the oscillations of cr^ are 
damped, due to the anharmonicities of the cosine 
potential. 
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FIG. 8. Eigenvalues X± (a) and o^/n (b). Nu- 
merical solutions for the variational ansatz (full line) 
and for the phase model (dotted line); N — 2 x 10 , 
g = 0.005 and r = 50. 

The conditions that preserve the superfluid state 
are that the final tunneling coupling JN is compa- 
rable to the on-site interaction and that the time 
scale of the process is slow enough to allow final 
small oscillations around the equilibrium. To de- 
stroy the phase coherence in this parameter range, 
one should raise the barrier faster in order to get 
larger oscilla tions, or raise it higher (case described 
in Sec.|VIA|). 



To summarise, in this Section we have deter- 
mined analytic expressions for the adiabaticity 
conditions for the internal dynamics, i.e. we iden- 
tified the time scale at which the barrier can 
be raised in order to obtain the same relative 
phase properties as expected for the ground state. 
Concerning a completed splitting process, a frag- 
mented condensate without any relative phase (no 
collapses and revivals) and characterised by a very 
narrow number distribution is reached in raising 
process with time scale r > 2r r . This means that 
the process has to be slower for large N and small 
g. Instead in an incomplete splitting process, the 
superfluid phase is preserved for r > 2ir/8JN[3. 
This time scale becomes shorter for large N. Both 
conclusions agree with the fact that in the Gross- 
Pitaevskii limit (N — > 00 and g — > 0) the conden- 
sate is phase coherent. 



VII. DISCUSSION AND CONCLUSIONS 

We have solved the two-mode model describing 
the splitting of a condensate by a potential barrier 
through a variational ansatz. We found coupling 
between the internal and the external dynamics of 
the mode functions. We have identified the regimes 
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in which the two dynamics decouple, and have con- 
cluded that in the case of splitting starting from a 
condensate in the ground state, they do not in- 
fluence each other in a dramatic way. Hence, the 
internal and external excitations created by raising 
the barrier can be estimated in a good approxima- 
tion independently and have been characterised as 
a function of the interaction strength, the number 
of atoms and the time scale of the process. 

From our analytic estimations, confirmed by nu- 
merical resuts, we were able to identify the time 
scales t z , t p and Tp, which define the adiabatic 
regime for the external and internal dynamics (re- 
spectively in the case of final fragmented or phase 
coherent condensate) 

T * >> h Tp>2Tr = 2 ik' T * > ^k- (38) 

It is interesting to compare the adiabaticity condi- 
tion for internal and external degress of freedom in 
the case of final fragmented condensate, when the 
splitting process can be consider to be completed. 
We normally found t z < r p ; the case N — 200, 
.9 = 5 sets a boundary where internal and external 
degrees of freedom enter simultaneously the adia- 
batic regime for r ~ 20 (see Fig. ||). To get t z > t p 
one needs gU\ to be large compared to the trap- 
ping frequency. The quantity gU\ is similar to the 
derivative of the chemical potential with respect 
to the total number of atoms d^/dN . From a 
Thomas-Fermi solution one gets that either in 1 
or 3 dimensions, it scales as a negative power of N 
and a positive power of g. So one needs to have 
very large g or very small N, and our model fails 
in both limits. Therefore, we claim that in the 
usual regime of many weakly interacting particles, 
the adiabaticity condition for the phase dynamics 
is more restrictive than the one for the external 
dynamics. 

We carried out a comparison of our model with 
the phase model, finding a substancial very good 
agreement. The numerical solution of the phase 
model consists in the integration of a time depen- 
dent Schrodinger equation. In practice, the num- 
ber of wavefunctions ip m (0) = exp(im0) / \/27r that 
one has to use increases linearly with N, which 
leads to numerical problems. In this sense for large 
N, the variational ansatz is more convenient and 
has proved to give reliable results. Moreover the 
variational ansatz treatment allows to include the 
external degrees of freedom in a natural way. 

In the previous results we did not take asymme- 
tries into account. It is possible to include them in 
our ansatz, through the unbalance in population 
mo and through non symmetrically centered wave 
functions. In the case of complete splitting, one 
ends up with a final constant unbalance in popu- 
lation. The phase coherence shows the only new 



feature that the center of the phase distribution is 
now drifting with a velocity f^i — /i2, where are 
the chemical potentials of the two separate conden- 
sates. A complete analysis of this case, which even 
consideres losses and fluctuations in the total num- 
ber of particles, can be found in |24|. Instead in the 
case of final phase coherent symmetric condensate, 
the asymmetry can destroy the phase coherence. 
The final unbalance in population might be so big, 
that the wavepacket describing the phase distribu- 
tion flies above the cosine potential: depending on 
the "kinetic energy" gUiuig, the cosine potential 
may become negligible and the same features (col- 
lapses, revivals) as in the fragmented condensate 
case are observed. 

Possible extensions of our model are the dynamic 
turning on of an optical lattice, where the initial 
harmonic trap is deformed into a many wells po- 
tential. The instantaneous version of such a pro- 
cess has recently been investigeted in 31 . Another 



problem of great significance is the inverse process, 
i.e. the merging of two condensates. This could al- 
low to refill a condensate and be an important step 
towards a continuous atom laser. 
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APPENDIX A: TWO-MODE 
VARIATIONAL ANSATZ APPROACH AND 
PHASE MODEL 

Dealing with two coupled condensates, in this 
paper we have often talked about the number dif- 
ference m and relative phase <fi. In this appendix, 
we will define number difference m and relative 
phase <j> operators, and derive the phase model 
Hamiltonian (^o|), discussing the approximations 
involved. 

We first define the operators 



1 



L x = 2 [ a i a 2 + a^ai J , 



Lz = 2 [ a 2 a 2 - a l a l J i 



(Ala) 
(Alb) 
(Ale) 
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such that the usual angular momentum commuta- 
tion relations are fulfilled [Li,Lj] = is^kL^.. After 
a small amount of algebra, the exact two-mode 
Hamiltonian in Eq.(|29|) can be rewritten as 

H = g{U! - 2U 2 )L 2 Z - 2JL X + gU 2 (L 2 X - L 2 ) . (A2) 

In the subspace of fixed even total number of atoms 
N, the spectrum of L z is given by all integer num- 
ber in the interval [—A/2, A/2]. The operator L z 
coincides with the number difference operator m. 
We have often treated the eigenvalues to as con- 
tinuous, but in general care should be exercised. 
Given the phase operator (j> such that [<f>, L z ] = i, 
it is well-known that in the phase-representation 
L z = rh = —id/dcf) and the eigenstates of L z (to) 
with eigenvalue to are (<fi\m) = exp(imcj)) / ^/2tt 
@- 

The state of the system |4>) can be in equivalent 
ways described as a superposition of eigenstates of 
to (c m is the number distribution) 



.|2l b i*)=x;^i( 

m 

N + 1 ^ 1 

2^ Cm o 



L+. +£_ \m 



(A7) 



2tt 



N/2 
m=-N/2 



(A3) 



or by a wave function in the (^-representation given 

by 

N/2 

= (<£|$) = —= c m exp(imcj}). (A4) 

V 27F m =-N/2 

Typical quantities characterizing |$) are 
- the width of the number distribution 



(m) = y^TO|c m | 

m 



u 2 = > m 2 \c m \ 2 



(A5) 



the width of the phase distribution 



(4>) = / 1<&(0)| 2 , 



4 = 



d<t>ci> 2 \®{4>)\ 2 -{4>f 



(A6) 



The uncertainty relations are the same one as for 
angular momenta |26| . 

In order to write the Hamiltonian in the phase- 
representation, we make some approximations 
which will lead to the simple phase model Hamilto- 
nian and discuss under which conditions such ap- 
proximations are valid. We will describe the pro- 
cedure only for the term L x , since for the term 
L 2 — L 2 an analogous one applies. Using the rais- 
ing and lowering operators L±, we write 



T-[^±i, , 



N + l 




1 - ( ^ 1 e"* 



V N+ 1 



We assume a narrow and centered number distri- 
bution (these assumptions will be quantified later 
on), so that we can expand the square roots at first 
order in (to/ N) 2 and get 



(4>\2L X \<5>) = 
= (N + l)cos 

+2(N+ 1) 



1/4 



(N + iy 



COS0 + 



m 



(N + iy 



■ sine 



-»b ~ Acos0$(0), 



(A8) 



where we can roughly estimate 

2 



^2 c m m 2 (4>\m) 



< 



(A9) 



< 



^ |c m |TO 2 TO" 



<^| Cm | 2 TO 2 ^ 

m m' 
(m)+er m 



to' 2 = 



0", , 



(to) 2 ) to 2 

ra— (m) — cr m 



(a m + (TO,) 2 ) O m < 



< (c m + I (to) I) a m , 
so that to/ 



< (cm + \ {fn)\) 2 y/c-m and where in 
a similar way \m\ = E m c m m(^|ra)| — (cr m + 
I (in) | ) \/f m . In an analogous way, the term L 2 —L 2 
gives 



4 ( L\ - I 2 



N 2 cos 20 -4 
N 2 cos 20. 



(A10) 



(m 2 + 1 ) cos 2(f) + 2m, sin 2(f> 



All together the neglected terms have to be smaller 
than all the other terms in the Hamiltonian, which 
implies 



Om + |(m)|) V^T« N 2 



2J - gNU 2 
N 



< g(U! - 2U 2 ). 



(All) 
(A12) 



For | (to) | <§; a m , which is the typical situation 
treated in this paper, from the first equation one 
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gets a m <C N 4 / 5 . In the opposite limit where 
| (to) | ^> a m instead one has 



I Ml 

N 



1/4' 



(A13) 



For a m — > this condition is not correct and be- 
comes simply | (m) | <C N. 

To summarise, in the specific case treated in 
this paper, where no unbalance of population be- 
tween the wells was assumed, the important condi- 
tion ofvalidity for the phase model Hamiltonian is 
Eq.( A12 ), which written for U2 — takes the form 



J < gNUx/2. 



(A14) 



We note that under such a condition, the ground 
state is characterised by a 4 > -y/l/iV or equiva- 
lently a m < y / WJl. 



APPENDIX B: COMPARISON BETWEEN 
VARIATIONAL ANSATZ AND PHASE 
MODEL 

In this subsection we will compare the phase 
model with our variational ansatz. For simplicity 
we set Ui = 0, but the following discussion could 
be repeated in the more general case. In particular 
for U2 = the phase model Hamiltonian reduces 
to the Josephson Hamiltonian pE 



d 2 

H = -gUi—r - JN cos< 



(Bl) 



The quantities gU± and JN are respectively the on- 
site energy splitting (charging energy) and the tun- 
neling coupling (Josephson coupling). It describe 
accurately the case of high barrier, where the two 
condensates are almost spatially separated, leading 
to a negligible U2, and are characterised by very 
small number fluctuations, making the Josephson 
Hamiltonian valid. The phase model Hamiltonian 



in Eq.(Bl) describes the motion of a particle in a 
cosine potential and the classical limit is obtained 
for cr m ,(j0 — > 0. In the following we discuss the 
classical and quantum limits comparing it directly 
with the variational ansatz results. 

In order to allow a c omp lete comparison with 
the Hamiltonian in Eq.(pTI), we now describe the 
coefficients c m in Eq.(^|) with four variational pa- 
rameters p, xi, toq and tp 



Af(p) exp 



~p +lXl 



too) 5 



(B2) 
exp [zto<p] , 



The variational parameters xi and tp are the vari- 
ables conjugate to p and mo, allow the dynamic 



evolution and vanish in the ground state. In the 
limit of broad number distribution this ansatz de- 
scribes a gaussian superposition of number states 
centered in to = too. The corresponding phase 
distribution $>[4>) is also a gaussian, centered in 
<fi = tp. Notice that in the case of symmetric split- 
ting treated before mo = and tp = Vi. 

Replacing again the sum over to from —N/2 to 
N/2 with an integral from —00 to 00 (for p > 1 
and xi such that < it), the widths of the num- 
ber and phase distribution are respectively given 
by Eqs.(|20|) and the expectation values (ajaj) and 
(a] 2 a 2 ) on the state 1$) are now 



+ N 
($K 2 ai, 2 |$) = — T m , (B3a) 

N 2 

($K? 2 a? >2 |*) = — + P TNm + m 2 , (B3b) 



($|4a 2 |$) - (*|4ai|*>* = 

N I a2 A <■ \ 

= y exp I — - I exp [tip) 



(B3c) 



1 2p m 2 ( \ 2 



a. Classical limit 



Let us first consider the classical limit and 
compare the results obtained with the variational 
ansatz with the phase model. We write the ex- 
pectation value of the Hamiltonian (Eq.(p9|)) for 
U2 = in the classical limit by setting a m = and 
0-0 = 



T~£cl = gUim — JNcostp; 



(B4) 



here we have assumed too <C N, neglected all terms 
o(l/N) <C 1 and all terms independent of too and 
p. One sees immediately that this exp ression co- 
incide with the Hamiltonian in Eq.(Bl), where the 
operators have been substituted with their expec- 
tation values ((to) = mo and (</>) = tp). So, the 
analogy with the phase model in the classical limit 
is straighforward po| . In particular if J > the 
stable equilibrium position is = 0; for a kinetic 
energy such that gU\m\ < JN the phase under- 
goes oscillations, otherwise if gU\m\ > JN the 
phase flies above the cosine potential. 



b. Quantum limit 

In the quantum limit the width of the number 
and phase distribution start to play an important 
role. The comparison between the two models now 
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is not so trivial, because we have on the one hand 
a wavefunction in the phase representation, $(0), 
and on the other hand 4 variational parameters 
(p, Xj, mo and tp) which follow a classical dynam- 
ics and should reproduce the features of $(</>). We 
consider here the simplified case of a wavefunction 
<&((f>) symmetric around = 0. In the variational 
ansatz picture, this corresponds to too = and <f = 
Vt. This last assumption is correct for symmetric 
initial conditions and preserved at all times, how 
can be checked esplicitly in the equations of mo- 
tion. 

Concerning the ground state of Hamiltonian 



(Bl), two different regimes can be identified: for 
gUi ^S> JN, the cosine potential can be neglected 
and the ground state is a flat wavefunction, which 
corresponds to a completely undefined relative 
phase between the two condensates. Instead for 
gUi -C JN the cosine potential is deep and the 
ground state is a localised wavefunction, which de- 
scribe a state with very well defined relative phase. 
In the variational ansatz approach, to find the 
ground state we set xi = 0, too = and ip = 0, and 
get H(p) = gU lP - JiVexp(-l/8p) [l - 2p/N 2 ]. 
For gUi <C JN, the minimum of this Hamiltonian 
is p — N/4 corresponding to the Gross-Pitaevkii 
limit ((70 ~ l/y/~N) and if the hopping decreases 
p — > (ct^ ^> 7r), reproducing the same results as 
the phase model. A quantitative comparison shows 
perfect agreement. 

The time evolution was discussed already in 
Sec. VI, where the analogy between the evolution 
of the phase distribution and the evolution 

of the variational parameters p and xi was car- 
ried out. The variational ansatz is able to pre- 
dict all the main features typical of the evolution 
of the phase distribution: collapses and revivals 
of the relative phase for gU\ ^> JN and breath- 
ing of the wavefunction in the cosine potential 
for JN ^S> gU\. Moreover the frequency of the 
small oscillations around the equilibrium point cal- 
culated in the variational approach concides with 
the splitting of the energy levels of the phase model 
Hamiltonian both in the limit of negligible cosine 
potential and in the limit of deep cosine potential. 



In this Appendix we have qualitatively com- 
pared the phase distribution described by the 
phase model with the Gaussian ansatz for the co- 
efficients of the state \<j>) in the number represen- 
tation. We have shown that in the classical limit 
the expectation values too and ip of number and 
phase (variational parameters indicating the cen- 
ters of the respective Gaussian distributions) are 
governed by the classical phase model Hamilto- 
nian. Moreover in the quantum description the 
parameters p and xi, related to the widths a m 



and tr (p , are able to reproduce the same features of 
the time evolution predicted by the phase model, 
whose limit of validity of the phase model (derived 
in App. 0) is J < gUxN/2. 

We discussed the limit of fragmented condensate 
(JN <§C gU\) and the limit of single phase coherent 
condensate (JN 3> gU\), As already mentioned, 
for fixed gN, the limit N — > oo corresponds to 
the Gross-Pitaevkii limit. Infact, unless J = 0, 
if N — > oo the condensate is always phase coher- 
ent. For finite N, the above discussion allows to 
judge whether the Gross-Pitaevskii description is 
still valid or not. 
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